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ABSTRACT 

We derive optimal periodic controls for entrainment of a 
self-driven oscillator to a desired frequency. The alternative ob- 
jectives of minimizing power and maximizing frequency range of 
lentrainment are considered. A state space representation of the 
oscillator is reduced to a linearized phase model, and the opti- 
mal periodic control is computed from the phase response curve 
using formal averaging and the calculus of variations. Compu- 
tational methods are used to calculate the periodic orbit and the 
phase response curve, and a numerical method for approximat- 
ing the optimal controls is introduced. Our method is applied 
to asymptotically control the period of spiking neural oscilla- 
tors modeled using the Hodgkin-Huxley equations. This example 
illustrates the optimality of entrainment controls derived using 
phase models when applied to the original state space system. 



1 INTRODUCTION 

The synchronization of oscillating systems is an important 
and extensively studied phenomenon in science, and also finds 
numerous engineering applications [IJ. Examples include the os- 
cillation of neurons ID, sleep cycles and other pacemakers in bi- 
ology mmiSl, semiconductor lasers in physics fS], and vibrating 
systems in mechanical engineering |7|. The asymptotic synchro- 
nization of an oscillator to a periodic control signal is called en- 
trainment, and is studied by examining the phase response curve 
(PRC) |[8][9l, which quantifies the shift in asymptotic phase due 
to an infinitesimal perturbation in the state. The classic phase 
coordinate transformation [ 10 1 for studying nonhnear oscillators 
was used together with formal averaging ifTTI to develop a model 
of coupled chemical oscillations ifTSl . Phase models are widely 
used in physics, chemistry, and biology fT3l to study systems 
where the phase, but not the state, can be observed, and where 
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the PRC can be approximated experimentally. Interest in con- 
trolling synchronization in electrochemical [14] and neural ifTSl 
systems has been increasing, and a method for approximating 
optimal waveforms for entrainment of phase-reduced oscillators 
by weak forcing has been proposed [16J . 



In this paper, we extend the theory of optimal entrainment of 
oscillators via weak, periodic controls llT6l to systems where the 
phase model has arbitrary PRC. We also present an efficient nu- 
merical method that accurately computes optimal waveforms by 
finding the maximum of a polynomial whose coefficients depend 
on the PRC of the entrained oscillator This enables an exami- 
nation of the important issue of how controls derived using the 
PRC perform when applied to entrain the associated oscillator in 
state space, which is the ultimate purpose of using phase models. 
In the following section, we discuss the phase coordinate trans- 
formation for a nonlinear oscillator and the available numerical 
methods for computing the PRC, and describe how averaging 
theory is used to study the asymptotic behavior of oscillating 
systems. In section 3, we use calculus of variations to derive 
theoretical entrainment controls that are optimal in the sense of 
minimum power or maximum entrainment range. The former 
is optimal when the natural frequency of the entrained oscillator 
is known to be either above or below the desired value, and the 
latter is useful when the natural frequency is in a neighborhood 
of the desired value, but unknown. We then present an efficient 
procedure for approximating these controls using Fourier series 
and Chebyshev polynomials. Finally in section 4, our approach 
is employed to entrain the Hodgkin-Huxley neuron model. The 
results suggest that optimal controls derived using a phase model 
are optimal for entrainment of the associated state space system. 



2 PHASE MODELS 

Consider a smooth ordinary differential equation system 



The weak ergodic theorem for measure-preserving dynamical 
systems on the torus [ 1 1 1 implies that for any (|), 



--f{x,u), x{Q)^xo, 



(1) 



where x{t) E R" is the state and u{t) e K is a control. Further- 
more, we require that ([T} has an attractive, non-constant limit 
cycle Y(f ) = Y(f + T), satisfying y = /(y, 0), on the periodic orbit 
r = {y € K" : y = y(0 for < f < T} C K". In order to study 
the behavior of this system, we reduce it to a scalar equation 



\}/= (£) + Z{\\t)u, 



(2) 



A(^) = (2(9 + ^)^(9)) 
1 r2n 
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exists as a smooth, 27t-periodic function in !f. By the formal 
averaging theorem [2\, the system 



which is called a phase model, where Z is the PRC and is 
the phase associated to the isochron on which x(f ) is located. The 
isochron is the manifold in R" on which all points have asymp- 
totic phase \|/(f) ifTTl . The conditions for validity and accuracy 
of this model have been determined 1 18 1, and the reduction is ac- 
complished through the well-studied process of phase coordinate 
transformation |fT9l . which is based on Floquet theory ll20]|2Tl . 
The model is assumed valid for inputs u{t) such that the solution 
x{t,xo,u) to ([T]i remains within a neighborhood of F. To com- 
pute the PRC, the period T —2k/(x) and the limit cycle y(r) must 
be computed to a high degree of accuracy. This is done using 
a method for determining the steady-state response of nonlinear 
oscillators ll22l based on perturbation theory ll23l and gradient 
optimization f24]. The PRC can then be computed by integrat- 
ing the adjoint of the linearization of ([T]l |25 1, or by using a more 
efficient and numerically stable spectral method developed more 
recently f26l. A software package called XPPAUT [27] is com- 
monly used by researchers to compute the PRC. We use a modi- 
fied spectral method in our implementation that is very accurate 
for stiff systems. 

Our goal is to entrain the system to a new frequency 
using a periodic control u{t) — k{0t) where k is 27t-periodic. We 
have adopted the weak forcing assumption, i.e. k ~ eki where 
ki has unit power, so the original system ([l} is guaranteed to 
traverse a neighborhood of F given this control. Now define a 
slow phase variable by (|)(f ) — \\i{t) — 0t, and call the difference 
Aco = CO — between the natural and forcing frequencies the fre- 
quency detuning. The dynamic equation for the slow phase is 



(j) = xj/ - = Aco + Z(0f + (j))fe(0f ) , 



(3) 



where (j) is called the phase drift. In order to study the asymptotic 
behavior of ^ it is necessary to eliminate the dependence on 
time, which can be accomplished by using formal averaging [ 12 1. 
Given a periodic forcing with frequency = 2k /T, we denote 
the forcing phase 9 = 0?. If ^ is the set of 27i-periodic functions 
on M, we can define an averaging operator (•) : M by 



1 /'2jr 

2-nL 



(4) 



(p = Aco + A((p) + O(e^) 



(6) 



approximates © in the sense that there exists a change of vari- 
ables (p = (|) + e/i((p,(|)) that maps solutions of Q to those of (|6]l. 
Therefore the weak forcing assumption k — eki with e << 1 al- 
lows us to approximate the phase drift equation by 



ip = Aco + A((p). 



(7) 



The averaged equation O is independent of time, and can be 
used to study the asymptotic behavior of the periodically forced 
system (|2]l where u = k{0t). 



3 ENTRAINMENT OF PHASE MODELS 

We call the system (|2]i entrained by a control u = k{0t) 
when the phase drift equation ^ satisfies cp = 0. This occurs 
when there exists a phase cp* satisfying Aco + A((p*) = 0, in 
which case the system is called entrainable. Defining the phases 
9- = argminip A((p) and (p+ = argmaX(p A((p), we can formulate 
entrainment as an optimal control problem. When the objective 
is to minimize the control power (^^), entrainability requires that 



Aco- 
Aco- 



^A((p+] 
'A((p-^ 



if 
if 



> CO, 
< CO. 



(8) 



We formulate the problem for > CO, and the case where < CO 
is symmetric. The constraint dH) can be added by adjoining it to 
the objective function using a multiplier X, resulting in 



mmj[k] = (yt^) -X(Aco + A((p+)) (9) 

= (k^)-X{A&+— Z(9 + (p+)/:(9)d9 
V 2k Jo 

= — / [fc(9)(/t(9)-A.Z(9 + (p+))-?iAco]d9 
2% Jo 



The Euler-Lagrange equation provides necessary conditions for 
the optimal solution, which is given by 



Because Z(9) is 27t-periodic, we represent it as a Fourier series, 



fe*(e) = -z(e + (p+). 

The constraint ([8]l can be used to solve for X, because 



1 r^^ X 

= Aco + A,(9+) =Aco+— / -Z(e + (p+)2de (10) 
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implies that X — — 2Aco/ (Z^). Consequently the minimum power 
control is 



Z(e) = xflo+ ^ fl„cos(ne) + b„sm{nQ), (14) 



^*(e) = -^z(e), 



(11) 



with power P — (Aco)^/ (Z^). We omit the phase ambiguity (p+ 
in the solution because entrainment is asymptotic. 

Now consider the dual problem where for fixed power P, 
a periodic waveform k{0t) is derived to maximize the locking 
range R[k] of natural frequencies CO for which the family of os- 
cillators {xj/ = co + Z(\(/)m : CO e (C0min,tt>max)} Can be enti-ained 
to a forcing frequency [16|. The locking range is given by 

R[k] = COmax - «min = ACOmin - ACOmax = A((p+) - A((p_ ), SO that 

adjoining the constraint on the power to the objective function 
using a multiplier X gives rise to the optimal control problem 

J[k;P]=R[k]-X{{k^)-P) (12) 
= A((p+)-A((p_)-M(fe2)_P) 

= (Z(e + (p+)fe(e)) - {Z{B + (?-)k{d))-X{{k^) -P) 

= / (fc(e)[z(e+(p+)-z(e + (p_)-U(e)]+A.p)de 

27t Jo 

Solving the Euler-Lagrange equation yields 

fe*(e)==i^[z(9+(p+)-z(e + (p_)]. 

The optimal solution k^, satisfies the constraint ^fe^) — f = 0, so 

1 /■27C / 1 \ 2 

2^/0 yrx) [z(6+^+)-z(6+'P-)]'de-i^ = o, 



and hence X=\ ^/qJp where Q = ( [Z(e + (p+) - Z(e + (p_)]2). 
Substituting this into (|5]l gives 

A*((p) = (z((p+e)^4e)) 
1 



2X 



(z((p + e)[z(e + (p+)-z(e + (p_)]) 



;i=l 



;i=l 



and we find that for (pi ,92 G [0,27t), 



(Z((pi + e)Z((p2 + e)> = \al + \Y^ {al + cos(«((pi - (pj)). 



Substituting this result into ( fTSt , we obtain 



(15) 



"^♦(^^ E(«» + ^n)['^os(n((p-(p+))-cos(n((p-(p_))]. 

(16) 

Let us denote the phase difference Acp = (p+ — (p_. Then 

e = ([z(0+(p+)-z(0+(p_)]2) 

= (z(e + 9+)2) - 2(z(e + (p+)z(e + (?-)> + (z(e + (p_)2) 

= EK' + ^')[l-cos(nA(p)]. 



(17) 



By substituting (p_ and (p+ into ( fTSI l, we obtain the optimal lock- 
ing range R[k^] as a function of A(p and the Fourier coefficients 
of Z, namely 

=A*((p+)-A*((p_) (18) 
= V^L(«^+^')[l-cos(nA(p)] = Vp^. 

n=\ 

Consequently, to find the optimal control k^, and the maximum 
locking range R\ki,\, it suffices to maximize Q in terms of A(p. 
The value of Acp that maximizes Q in ( [TtI i also satisfies the first 
order condition 2'(A(p) = Lir=i "("^n + ^«) sin(nA(p) = 0, hence 
there exists a "generic" solution Acp = 7t, which may not be opti- 
mal. Observe that if we set y = cos(Acp), then 

e(Acp) = q{y) = f {al + bl)[l - r„(y)], (19) 

where T„ is the Chebyshev polynomial of the first kind. 
Therefore a straightforward criterion for the existence of superior 
solutions is to check whether q attains its supremum on (—1,1). 
In that case we choose Acp — ± arccos(3'*), and otherwise we 
choose Acp = K. The optimal waveform is given by 



kM = VP/Q[Z{Q + A^,) -z(e)]. 



(20) 



Vp7e(z(cp + e)[z(e + cp+)-z(e + cp_)]). (i3) 



We omit the phase ambiguity cp_ in (|20] | because entrainment is 
asymptotic. The two possible values for A(|) result in two optimal 
solutions when the criterion for (IT9^ holds. 



4 ENTRAINMENT OF NEURONS 

The notion of modeling the dynamics of neurons in the hu- 
man brain as oscillators has gained wide acceptance among re- 
searchers in neuroscience and mathematical biology [ 17, 28 1. Be- 
cause the ability to control the synchronization of neural dynam- 
ics has important research and clinical implications II29I I30L it is 
important to explore the pertinence of the entrainment paradigm 
to neural systems. We consider the entrainment of a neuron by an 
external stimulus, and use as an example the model of Hodgkin 
and Huxley [31 1. Starting with the commonly used parameteri- 
zation flTl, we reduce the system to the phase model and com- 
pute optimal entrainment controls. The objective is either to en- 
train the model to a given frequency with minimum power (|9), 
or to maximize the range of frequencies (and hence the number 
of neurons) that can be entrained by a control of fixed power 
(fT2l i. For a given waveform k{0t) where is in a neighborhood 
of the natural frequency co, we can numerically approximate the 
power actually required for entrainment. This allows us to com- 
pute the approximately triangular region of entrainability called 
the Arnold tongue, which is the plot of the minimum amplitude 
^/P required for entrainment versus forcing frequency 0, and 
which is commonly used to visualize the asymptotic properties 
of an oscillating system 1.13. .32J . This will be used to illustrate 
the performance of the controls that we have derived. 

The Hodgkin-Huxley model describes the propagation of ac- 
tion potentials in neurons, specifically the squid giant axon, and 
is used as a canonical example of neural oscillator dynamics. The 
equations are 



cV 


= h + I{t)- 


- gj,My - ^NaW - gK{y - - gL(y 


m 


= a,n{V){l- 


- m) - bm{V)m, 


h 


= ah{V){\- 


-h)-bh{V)h, 


n 


= fl„(y)(l- 


-n)-b„{V)n, 




am{y) = 


0.1(y+40)/(l-exp(-(y+40)/10)), 






4exp(-(y + 65)/18), 




ah{y) = 


0.07exp(-(y + 65)/20), 




bh{V) = 


l/(l+exp(-(y + 35)/10)), 




a„{V) = 


0.01(y + 55)/(l-exp(-(y + 55)/10)), 




bn{y) = 


0.125exp(-(y + 65)/80). 



(21) 

The variable V is the voltage across the axon membrane, and 
m, h, and n are the ion gating variables. is a baseline cur- 
rent that induces the oscillation, and I{t) is the control input. 
The units of V are millivolts and the units of time are millisec- 
onds. We analyze this system of differential equations as an os- 
cillator i — f{x,u), with a periodic limit cycle y(f) = y(f + T) 
present when u = 0. Using the standard parameters Vnu = 
50 mV, Vk = -77 mW, Vl = -54.4 mV, = 120mS/cm2, 
= 36 mS/cm^, = 0.3 mS/cm^, 4 = lO^A/cm^, and c = 
1 juF/cm^, we compute the limit cycle, which is shown for 
the voltage V in Figure [T] The period is computed as T = 
14.63842 ± 10^^ ms. The "spiking" behavior of the oscillator 
indicates that this system is stiff, and hence ill-conditioned for 
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Figure 1 . Hodgekin-Huxley limit cycle (left) and "spiking" (right) 
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Figure 2. Hodgekin-Huxley phase response curve (PRC) 

numerical integration. We use a second order Adams-Bashforth 
solver to integrate these equations with a relative error tolerance 
of 10^^. The PRC is computed along the limit cycle with an ini- 
tial conditionxo = {Vo,mQ,no,ho) = (0,0.51916,0.2999,0.4812) 
corresponding to V(/(0), and the result is shown in Figure |2] An 
absolute error lower than 10^^ is maintained by using a grid with 
step size 0.002. The first and second zero crossings occur at 
\(/ ~ 0.4617 and \\f = 4.2242, respectively. Note that u is least 
effective at the start of the cycle, when the neuron is spiking. 

We expand the PRC in a Fourier series as 
in (fT4l i by using the discrete Fourier transform of 
{Z{j) : j ^2% j/N, ;■ = !,... 5000} to approximate the 
coefficients. This gives us Z(n) = ^^^j Z(j)a)^ ''where 

= e^^Jii/A'^ tjjg estimates are a„ ^ 5R(Z(n)) • 2/N and 
b„ = — 3(Z(«)) - l/N. Because of the phase ambiguity, the 
choice of XQ e F that is used to compute y(f) influences the 
values of a„ and bn, but not the value of \a„ + ibn\. We take 20 
Fourier modes for our approximation. The total power of the 
Hodgkin-Huxley PRC as a periodic waveform is 0.0387, and the 
modes k=\,2,...,5 have power 0.01706, 0.01649, 0.00473, 
0.00048, and 0.00001, respectively. The modes 2 and 3 have 
significant power, hence it is insufficient to use a single mode 
to approximate the PRC. The minimum power waveform (fTTt 
is a re-scaled PRC. To compute the maximum range waveform 
( l20b . we find that a value of y^, = —0.05287 maximizes the 
polynomial q in ( fT9] ) on (—1,1), hence the "generic" solution 
A(p = 7t is not optimal, so we use Acp ~ arccos(x*) = 1.623690 
and get Q w 0.10976 ± 1. The polynomial q, its maximum, and 
the maximum range control waveform (l20t with unity power are 
shown in Figure [3] 




Figure 3. q{x) (left) with A(p (marl<ed), and max range k^. 




Figure 4. Arnold tongues for Hodgkin-Huxley phase model |2): Minimum 
power theory (dashed line) and computation (o); Maximum range theory 
(solid line) and computation (•). The minimum power control functions 
as intended only to increase frequency while the maximum range control 
has a useful symmetry property. 

To evaluate the entrainability of a phase-reduced system by 
a given waveform, we compute the Arnold tongue by deter- 
mining the power required for entrainment at a given frequency 
0. The key idea is that if entrainment does indeed occur, then 
the response of the oscillator is periodic with a period equal to 
T — iTi/Q. If the solution to (|2]i with u ~ k{0t) is sampled at this 
interval and the sequence {v(y7")};GN converges, it follows that 
the control u entrains the phase model. We determine the power 
f * (0) required for the sequence to converge by performing a bi- 
section search, using 150 points of the sequence as a test. A plot 
of ^/l\{0) vs. generates the resulting Arnold tongue. The 
distinction between the solutions ( fTTT i and (|20] | obtained by us- 
ing the alternative objectives is illustrated in Figure^] The results 
for (fTTT l on the irrelevant range are omitted in other figures. The 
Arnold tongues for the phase reduced system are presented in 
Figure |5] Note that the actual Arnold tongues are not linear, and 
the required power to decrease (increase) the frequency is lower 
(higher) than predicted by the theory. An issue of fundamental 
importance is how well the entrainment control works when it is 
applied to the original Hodgkin-Huxley system. Figure |6] shows 
P* (0) VS. when the same control waveforms are applied to 
the original system i2T[ . The power required to entrain the state 
space model to a frequency CO is similar to the theoretical pre- 
diction near the natural frequency. By comparing Figures 2]and 
|6] one sees that the relative entrainability of the phase and state 
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Figure 5. Arnold tongue for Hodgkin-Huxley phase model j2): Mini- 
mum power theory (dashed line) and computation for increase (o) and 
decrease (*) of frequency; Maximum range theory (solid line) and com- 
putation (•); sine wave computation (+). The minimum power waveform 
for increasing (o) (decreasing {*)) CO matches the theory (dashed line) 
closely near COq for CO > COq (CO < COq). Similarly, the maximum range 
waveform (•) matches the theory (solid line) closely near cOo, and can be 
effectively applied to increase or decrease the frequency The sine wave 
(+) has the worst performance. 
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Figure 6. Arnold tongue for Hodgkin-Huxley state-space model (21]: 
Computed minimum power control for increase (o) and decrease (*) of 
frequency, and theory (dashed line); Computed maximum range (•), and 
theory (solid line); sine wave (+). 

models by the tested waveforms is nearly identical for values of 
near the natural frequency cOq. This is strong evidence that op- 
timal entrainment waveforms for a phase-reduced oscillator (|2]l 
are optimal in the same sense for the state-space system ([T) from 
which the reduced model is derived. 



Conclusions 

We have presented a method for optimal entrainment of 
oscillators given the alternative objectives of minimum control 
power and maximum range of entrainability. The method that 
we derived is based on the phase response curve of the oscilla- 
tor and formal averaging theory. We examine the entrainment 



of phase-reduced Hodgkin-Huxley neurons as an example prob- 
lem, and compute Arnold tongues to evaluate the effectiveness of 
our controls. Their performance closely matches the theoretical 
bounds when the weak forcing requirement is fulfilled. The opti- 
mal waveforms produce a similar result when applied to the orig- 
inal model, which suggests that optimal entrainment controls for 
a phase model are optimal for the original system, provided the 
oscillator remains within a neighborhood of its limit cycle. This 
work provides a basis for evaluating the effectiveness of phase 
reduction techniques for the control of oscillating systems. The 
approach described is of direct interest to researchers in chem- 
istry and neuroscience, and may also be applied to vibration con- 
trol in engineered systems. 
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